Explicit Solution is available

source("../R/exppop.R")

exppop

## function (T, P0, r, K) 
## {
##     P = P0 * exp(r * T)
##     if (P > K) {
##         P = K
##     }
##     return(P)
## }

# gives population after any time given an initial population

# 20 rabbits, growth rate of 0.01 how many in 30 years
exppop(T=30, P0=20, r=0.01, K=1000)

## [1] 26.99718

# if we want to see how population evolves over time - generate a time series by running our model for each point in time

initialrabbits = 20
years = seq(from=1, to=100, by=2)
Ptime = years %>% map_dbl(~exppop( P0=initialrabbits, r=0.01, K=1000, T=.x))

# keep track of what times we ran
Ptime = data.frame(P=Ptime, years=years)

ggplot(Ptime, aes(years,P))+geom_point()+labs(x="years",y="Rabbit Population")

Difference Equations

Population models can be discrete (rather than continuous)

So we could implement them as difference equations and iterate

source("../R/discrete_logistic_pop.R")
# notice how a for loop is used to iterate

# how many rabbits after 50 years given a growth of 0.1
# starting with 1 rabbit - but a carrying capcity of 500

discrete_logistic_pop

## function (P0, r, K, T = 10) 
## {
##     pop = P0
##     for (i in 1:T) {
##         pop = pop + r * pop
##         pop = ifelse(pop > K, K, pop)
##     }
##     return(pop)
## }

discrete_logistic_pop(P0=1, r=0.05, K=200, T=50)

## [1] 11.4674

# save results
discrete_result = discrete_logistic_pop(P0=1, r=0.05, K=200, T=50)

# lets also keep the parameters for use later
P0=1
r=0.05
K=200
T=50

Compare discrete and analytic results

source("../R/exppop.R")

exppop(P0=P0, r=r, K=K, T=T)

## [1] 12.18249

analytic_result = exppop(P0=P0, r=r, K=K, T=T)

analytic_result

## [1] 12.18249

discrete_result

## [1] 11.4674

# why are they different
# look at trajectories

growth_result = data.frame(time=seq(from=1,to=100))

growth_result$Panalytic = growth_result$time %>% map_dbl(~exppop( P0=1,r=0.05, K=200,T=.x ))

growth_result$Pdiscrete = growth_result$time %>% map_dbl(~discrete_logistic_pop( P0=1,r=0.05, K=200,T=.x ))

tmp = growth_result %>% gather(key="Ptype",value="P",-time)
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()

# try running them for longer time periods to see what happens 
# change the value of r, K , P0 - see how it effects the results

Solving using numeric integration

Using a solver….when you can’t do the integration by hand :)

For example, if you added a carrying capacity as threshold where it stops growing

In this case we integrate by iteration - approximates analytic integration

Numerical integration with ODE

Implement the differential equation as a function that

returns the derivative (as a list)

*inputs time, the variable(s) and a parameter list

(it needs time even though you don’t use it)

My convention: name derivative functions starting with d to remind myself that they are computing a derivative

ODE

Only works for Ordinary Differential Equations - single independent variable (in our case time)

Partial differential equations - more than 1 independent variable (e.g x and y if changing in space)

R has a solver called ODE for solving ordinary differential equations from package desolve

ODE requires

initial conditions
values of independent where we want values of dependent variable (e.g times where we want population)
the derivative as a function
a list that contains all parameter values (or if you have only one parameter then you can use a single value)

source("../R/dexppop.R")

dexppop

## function (time, P, r) 
## {
##     dexpop = r * P
##     return(list(dexpop))
## }

library(deSolve)
initialrabbits = 20
years = seq(from=1, to=100, by=2)

# run the solver
Ptime = ode(y=initialrabbits, times=years, func=dexppop,parms=c(0.01))
head(Ptime)

##      time        1
## [1,]    1 20.00000
## [2,]    3 20.40404
## [3,]    5 20.81623
## [4,]    7 21.23675
## [5,]    9 21.66576
## [6,]   11 22.10344

colnames(Ptime)=c("year","P")

# notice that there are additional pieces of information year, including the method used for integration
attributes(Ptime)

## $dim
## [1] 50  2
## 
## $dimnames
## $dimnames[[1]]
## NULL
## 
## $dimnames[[2]]
## [1] "year" "P"   
## 
## 
## $istate
##  [1]   2  52 105  NA   6   6   0  36  21  NA  NA  NA  NA   0   1   1  NA  NA  NA
## [20]  NA  NA
## 
## $rstate
## [1]  2.00000  2.00000 99.98839  0.00000  1.00000
## 
## $lengthvar
## [1] 1
## 
## $class
## [1] "deSolve" "matrix" 
## 
## $type
## [1] "lsoda"

# this also means you need to extract just the data frame for plotting
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")

# this also works (of course function can be by order)
Ptime=ode(initialrabbits, years, dexppop,0.03)
colnames(Ptime)=c("year","P")
ggplot(as.data.frame(Ptime), aes(year,P))+geom_point()+labs(y="Population", "years")

additional parameters

You can play a bit with changing your function to something that you can’t integrate “by hand”

BUT we might want more parameters

to work with ODE, parameters must all be input as a single list; similar to how we return multiple outputs from a function

see example below..lets add a carrying capacity

R code with carrying capacity

source("../R/dexppop_play.R")

dexppop_play

## function (time, P, parms) 
## {
##     dexpop = parms$r * P
##     dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
##     return(list(dexpop))
## }

# create parameter list
initalrabbits=2
newparms = list(r=0.03, carry_capacity=300)

#apply solver
results=ode(initialrabbits, years, dexppop_play,newparms)
head(results)

##      time        1
## [1,]    1 20.00000
## [2,]    3 21.23675
## [3,]    5 22.54993
## [4,]    7 23.94435
## [5,]    9 25.42499
## [6,]   11 26.99719

# add more meaningful names
colnames(results)=c("year","P")

# plot
ggplot(as.data.frame(results), aes(year,P))+geom_point()+labs(y="Population", "years")

# try again with different parameters
alternativeparms = list(r=0.04, carry_capacity=500)
results2=ode(initialrabbits, years, dexppop_play,alternativeparms)

## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 81.5108
##

## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps

## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go

# look at results
head(results2)

##      time        1
## [1,]    1 20.00000
## [2,]    3 21.66575
## [3,]    5 23.47022
## [4,]    7 25.42499
## [5,]    9 27.54257
## [6,]   11 29.83650

colnames(results2)=c("year","P_parm2")

# plot
ggplot(as.data.frame(results2), aes(year,P_parm2))+geom_point()+labs(y="Population", "years")

# compare by combining into a single data frame
both = inner_join(as.data.frame(results), as.data.frame(results2))

## Joining with `by = join_by(year)`

both_p = both %>% gather(key="model", value="P", -year)
ggplot(both_p, aes(year,P, col=model))+geom_point()+labs(y="Population", "years")

# try playing on your own - modify the function in some way

Differential Equation, Difference (Iteration by hand) comparison

Remember we have 3 ways now to calculate population

analytical solution - based on integration (exppop.R) BEST

using an ode solver for numerical approximation (exppop_play.R)

numerical integration using in discrete steps (discrete_logistic_pop.R)

Compare analytical, difference and ODE

Finally look at continuous derivative using ODE solve Needs initial conditions differential equation *parameters

source("../R/dexppop_play.R")

dexppop_play

## function (time, P, parms) 
## {
##     dexpop = parms$r * P
##     dexpop = ifelse(P > parms$carry_capacity, 0, dexpop)
##     return(list(dexpop))
## }

# set up using the same parameters
pcompare = list(r=r,carry_capacity=K)


# now run our ODE solver
result = ode(P0, growth_result$time, dexppop_play, pcompare)
head(result)

##      time        1
## [1,]    1 1.000000
## [2,]    2 1.051273
## [3,]    3 1.105172
## [4,]    4 1.161835
## [5,]    5 1.221404
## [6,]    6 1.284027

# we already have time - so just extract population
growth_result$Pdifferential = result[,2]

# compare all 3 approaches
tmp = growth_result %>% pivot_longer(cols=-time,names_to="Ptype",values_to="P")
ggplot(tmp, aes(time,P, col=Ptype))+geom_point()

# notice Pdifferential is closer to Panalytic than Pdiscrete

Other Examples

All differential and difference equations are approximations The analytical solution is exact

Notice that differential equations is a bit more accurate!

Diffusion

Diffusion can be implemented as a partial differential equation

Complicated to solve - but there are tool in R for specific types of partial differential equations

More info on differential equations in R

Diffusionn would require partial derivatives - time and space! it gets much more tricky …beyond this course

But we can appoximate diffusion with a difference equation - and iterative to get an estimate of how diffusion works

Example of Diffusion - difference equation implementation to see what some issues can be

Diffusion in 1 dimension

Diffusion in one dimension through time

Diffusion data structure

R implementation

source("../R/diffusion.R")

# run our diffusion model (iterative difference equation) with initial concentration of 10, for 8 timestep (size 1m), and 10 space steps (size 1s)
# using diffusion parameters 0.5 s/m2, 10 m2
result = diff1(initialC=10, nx=10, dx=1, nt=8, dt=1, D=0.5, area=10)

# a list is returned with our 3 data frames for concentration (conc), qin and qout
result

## $conc
##           [,1]     [,2]     [,3]      [,4]      [,5]        [,6]        [,7]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [2,]  7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000 0.000000000
## [3,]  6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000 0.000000000
## [4,]  5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000 0.000000000
## [5,]  4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000 0.000000000
## [6,]  4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625 0.000000000
## [7,]  4.189453 3.142090 1.745605 0.6982422 0.1904297 0.031738281 0.002441406
## [8,]  3.927612 3.054810 1.832886 0.8331299 0.2777100 0.064086914 0.009155273
##              [,8] [,9] [,10]
## [1,] 0.0000000000    0     0
## [2,] 0.0000000000    0     0
## [3,] 0.0000000000    0     0
## [4,] 0.0000000000    0     0
## [5,] 0.0000000000    0     0
## [6,] 0.0000000000    0     0
## [7,] 0.0000000000    0     0
## [8,] 0.0006103516    0     0
## 
## $qout
##           [,1]     [,2]     [,3]     [,4]       [,5]       [,6]        [,7]
## [1,] 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [2,] 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
## [3,]  7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000 0.000000000
## [4,]  5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000 0.000000000
## [5,]  4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000 0.000000000
## [6,]  3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406 0.000000000
## [7,]  2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219 0.006103516
## [8,]  0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000 0.000000000
##      [,8] [,9] [,10]
## [1,]    0    0     0
## [2,]    0    0     0
## [3,]    0    0     0
## [4,]    0    0     0
## [5,]    0    0     0
## [6,]    0    0     0
## [7,]    0    0     0
## [8,]    0    0     0
## 
## $qin
##      [,1]      [,2]     [,3]     [,4]     [,5]       [,6]       [,7]
## [1,]    0 25.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
## [2,]    0 12.500000 6.250000 0.000000 0.000000 0.00000000 0.00000000
## [3,]    0  7.812500 6.250000 1.562500 0.000000 0.00000000 0.00000000
## [4,]    0  5.468750 5.468750 2.343750 0.390625 0.00000000 0.00000000
## [5,]    0  4.101562 4.687500 2.636719 0.781250 0.09765625 0.00000000
## [6,]    0  3.222656 4.028320 2.685547 1.074219 0.24414062 0.02441406
## [7,]    0  2.618408 3.491211 2.618408 1.269531 0.39672852 0.07324219
## [8,]    0  0.000000 0.000000 0.000000 0.000000 0.00000000 0.00000000
##             [,8] [,9] [,10]
## [1,] 0.000000000    0     0
## [2,] 0.000000000    0     0
## [3,] 0.000000000    0     0
## [4,] 0.000000000    0     0
## [5,] 0.000000000    0     0
## [6,] 0.000000000    0     0
## [7,] 0.006103516    0     0
## [8,] 0.000000000    0     0

# used filled contour to plot results
head(result$conc)

##           [,1]     [,2]     [,3]      [,4]      [,5]        [,6] [,7] [,8] [,9]
## [1,] 10.000000 0.000000 0.000000 0.0000000 0.0000000 0.000000000    0    0    0
## [2,]  7.500000 2.500000 0.000000 0.0000000 0.0000000 0.000000000    0    0    0
## [3,]  6.250000 3.125000 0.625000 0.0000000 0.0000000 0.000000000    0    0    0
## [4,]  5.468750 3.281250 1.093750 0.1562500 0.0000000 0.000000000    0    0    0
## [5,]  4.921875 3.281250 1.406250 0.3515625 0.0390625 0.000000000    0    0    0
## [6,]  4.511719 3.222656 1.611328 0.5371094 0.1074219 0.009765625    0    0    0
##      [,10]
## [1,]     0
## [2,]     0
## [3,]     0
## [4,]     0
## [5,]     0
## [6,]     0

filled.contour(result$conc, xlab="Time", ylab="Distance")

# or if you prefer this orientation (Distance on x axis)
filled.contour(t(result$conc), ylab="Time", xlab="Distance")

Change parameters (diffusivity D, and space and time steps (dx, dt))

# changes diffusivity and other parameters particularly
# diffusivity, dx and dt
res=diff1(initialC=10,nx=10,dx=1,nt=10,dt=30,D=0.006,area=1)

filled.contour(res$conc, xlab="Time", ylab="Distance")

# we can also see how much material moved from place to place each time step
filled.contour(res$qin, xlab="Time", ylab="Distance")

# play with time step, space step and parameters

Play

Try running the diffusion model with different time steps, space steps and parameters

Some examples with different parameters and space/time steps

# what if we increase diffusivity
resfast=diff1(initialC=10,nx=10,dx=0.5,nt=10,dt=10,D=0.08,area=1)
filled.contour(resfast$conc, xlab="Time", ylab="Distance")

# Discretization Issue Example
resunstable=diff1(initialC=10,nx=10,dx=1,nt=10,dt=10,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="Time (fraction of hour)",ylab="Distance Along Path (m)", main="Pollutant Diffusion")

# this illustrates the problem with difference equations (and the challenges that methods for numerical integration try to overcome)
# if things are changing quickly we need to use much smaller time, space steps to avoid overshoot and instability

# so lets cut our step size by 10 (dt) (but then add 10 more steps (nx) to cover the same distance)
resunstable=diff1(initialC=10,nx=100,dx=1,nt=10,dt=1,D=0.8,area=1)
filled.contour(resunstable$conc, xlab="time",ylab="Distance Along Path")

More on Dynamic Models

Integration